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ABSTRACT 

The  Lagrangian  predictability  of  general  circulation  models  is  limited  by  the  need  for  high-resolution  data 
streams  to  constrain  small-scale  dynamical  features.  Here  velocity  observations  from  Lagrangian  drifters 
deployed  in  the  Gulf  of  Mexico  during  the  summer  2012  Grand  Lagrangian  Deployment  (GLAD)  experiment 
are  assimilated  into  the  Naval  Coastal  Ocean  Model  (NCOM)  4D  variational  (4DVAR)  analysis  system  to 
examine  their  impact  on  Lagrangian  predictability.  NCOM-4DVAR  is  a  weak-constraint  assimilation  system 
using  the  indirect  representer  method.  Velocities  derived  from  drifter  trajectories,  as  well  as  satellite  and  in 
situ  observations,  are  assimilated.  Lagrangian  forecast  skill  is  assessed  using  separation  distance  and  angular 
differences  between  simulated  and  observed  trajectory  positions.  Results  show  that  assimilating  drifter 
velocities  substantially  improves  the  model  forecast  shape  and  position  of  a  Loop  Current  ring.  These  gains  in 
mesoscale  Eulerian  forecast  skill  also  improve  Lagrangian  forecasts,  reducing  the  growth  rate  of  separation 
distances  between  observed  and  simulated  drifters  by  approximately  7.3  km  day-1  on  average,  when  com¬ 
pared  with  forecasts  that  assimilate  only  temperature  and  salinity  observations.  Trajectory  angular  differ¬ 
ences  are  also  reduced. 


1.  Introduction 

The  importance  of  Lagrangian  forecasts  was  seen 
in  the  wake  of  the  Deepwater  Horizon  oil  spill  in  the 
northern  Gulf  of  Mexico  (GOM).  The  agencies 
responding  to  this  disaster  required  information  about 
the  transport  and  dispersion  of  water  parcels  and  in  turn 
oil.  Unfortunately,  while  numerical  models  and  data 
assimilation  methods  continue  to  improve,  trajectory 
forecasts  remain  challenging. 

A  number  of  studies  have  investigated  Lagrangian 
predictability  using  model  velocities  and  observed 
drifter  trajectories  (Thompson  et  al.  2003;  Ozgokmen 
et  al.  2003;  Barron  et  al.  2007;  Huntley  et  al.  2011).  These 
studies  show  that  model  trajectory  uncertainties  typi¬ 
cally  grow  rapidly,  often  on  the  order  of  0.5-1  km  h_1. 
This  is  likely  due  to  small  model  velocity  uncertainties 
that  accumulate  during  trajectory  integration  (Huntley 
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et  al.  2011).  Model  forecast  velocities  remain  prone  to 
errors  due  to  the  lack  of  observations  available  for  as¬ 
similation  and  validation.  Increasing  model  resolution 
aggravates  this  problem,  since  even  more  observations 
are  needed  to  constrain  the  smaller  scales  that  are  re¬ 
solved.  As  more  observations  become  available  in  the 
future,  the  constraints  they  provide  through  assimilation 
should  directly  translate  into  improved  trajectory  fore¬ 
cast  skill. 

The  types  of  observations  assimilated  and  their  spatial 
and  temporal  resolution  also  play  an  important  role. 
Temperature,  salinity,  and  sea  surface  height  (SSH, 
measured  along-track  by  satellite  altimeters)  observa¬ 
tions  are  typically  assimilated  in  operational  ocean 
models.  However,  it  should  be  noted  that  along-track 
SSH  observations  are  not  assimilated  directly;  they  are 
used  to  calculate  synthetic  profiles  of  temperature  and 
salinity,  and  it  is  these  synthetic  profiles  that  are  assim¬ 
ilated.  Since  these  observations  directly  constrain  only 
the  model  mass  field,  thus  influencing  only  the  baro- 
tropic  component  of  the  flow,  their  impact  on  trajectory 
forecast  skill  is  unclear.  In  addition,  SSH  measurements 
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are  highly  resolved  in  space  and  time  along  the  satellite 
track,  but  are  coarsely  resolved  in  the  cross-track  di¬ 
rection  with  track  repeat  intervals  as  long  as  10  days. 
When  assimilated,  the  disparity  in  space  and  time  scales 
of  the  altimeter  measurements  lead  to  biases  in  the 
mesoscale  forecasts  that  are  very  difficult  to  quantify 
(Jacobs  et  al.  2014a).  These  biases  certainly  impact  the 
skill  of  model  trajectory  forecasts. 

When  adding  assimilation  of  Lagrangian  trajectories 
(or  their  derived  velocities)  to  the  typical  assimilation 
constraints  on  the  model  mass  fields,  the  trajectory 
forecast  skill  improves.  A  variety  of  methods  have  been 
explored  for  assimilating  Lagrangian  information  into 
ocean  models.  Toner  et  al.  (2001)  used  a  normal  mode 
analysis  method  with  constrained  optimization  to  mini¬ 
mize  changes  in  model  velocities  while  matching  ve¬ 
locities  along  trajectories.  Kuznetsov  et  al.  (2003)  and 
Ide  et  al.  (2002)  presented  a  method  for  the  direct  as¬ 
similation  of  Lagrangian  data  using  augmented  tracer 
advection  equations  that  track  the  correlations  between 
the  flow  and  the  tracers  using  an  extended  Kalman  filter. 
Molcard  et  al.  (2005)  used  a  statistical  method  to  cor¬ 
relate  model  and  drifter  velocities.  Taillandier  et  al. 
(2006)  describe  the  Lagrangian  variational  analysis 
(LAVA)  method  that  minimizes  the  distance  between 
observed  and  model  drifter  positions.  Additionally, 
Nodet  (2006)  presents  an  identical  twin  experiment  us¬ 
ing  four-dimensional  variational  data  assimilation 
(4DVAR)  in  an  idealized  wind-driven  midlatitude  box 
model.  These  previous  studies  focused  on  sequential 
methods  or  variational  techniques  applied  to  simplified 
or  idealized  ocean  models  with  simulated  data. 

Using  synthetic  data  from  a  model  control  run, 
Vernieres  et  al.  (2011)  showed  that  the  assimilation  of 
trajectory  data  is  surprisingly  effective  in  constraining 
the  model  SSH  field  to  capture  eddy  shedding  events, 
much  more  so  than  the  assimilation  of  comparable 
Eulerian  velocity  data.  Assimilating  trajectory  positions 
carries  with  it  the  challenge  of  estimating  appropriate 
covariances.  This  problem  is  circumvented  if  derived 
velocities  are  used  instead.  Carrier  et  al.  (2014)  dem¬ 
onstrated  significant  improvements  in  forecast  skill  in 
model  velocity,  salinity,  and  SSH  fields  when  velocities 
derived  from  drifters  were  assimilated  in  addition  to 
Eulerian  observations.  This  study  goes  a  step  further  and 
seeks  to  quantify  how  Lagrangian  forecast  skill  is  af¬ 
fected  by  the  assimilation  of  drifter-derived  velocities. 
For  this  purpose,  the  series  of  experiments  conducted  by 
Carrier  et  al.  (2014),  using  a  dynamically  consistent  4D 
variational  approach  to  assimilate  velocities  derived 
from  roughly  300  observed  drifters  in  the  northern 
GOM  into  an  operational  ocean  model,  is  analyzed  for 
their  skill  in  predicting  future  drifter  trajectories. 


Additionally,  a  first  step  toward  determining  the  spatial 
extent  of  these  Lagrangian  prediction  improvements  is 
undertaken  with  a  drastic  data-denial  experiment. 

This  paper  is  organized  as  follows.  In  section  2 
the  model  and  assimilation  system  are  described.  In 
section  3  the  observations  and  the  data  processing  are 
discussed.  The  methods  used  computing  model  trajec¬ 
tories  and  their  comparison  with  observations  is  dis¬ 
cussed  in  section  4.  Trajectory  forecast  skill  results  are 
described  in  section  5.  Finally,  section  6  provides  a 
summary  and  conclusions. 

2.  Model  and  data  assimilation  system 

The  model  used  here  is  a  GOM  regional  version  of 
the  U.S.  Navy  Coastal  Ocean  Model  (NCOM),  a  free- 
surface  primitive  equation  model  that  uses  hybrid  gen¬ 
eralized  vertical  coordinates  (Martin  2000).  NCOM  is 
used  operationally  to  produce  forecasts  of  temperature, 
salinity,  velocity,  and  SSH  for  numerous  global  and  re¬ 
gional  applications. 

The  GOM  regional  NCOM  domain  spans  18°-31°N 
latitude  and  79°-98°W  longitude  with  a  horizontal 
resolution  of  6  km.  Of  the  50  vertical  layers,  the  top  10 
are  free  sigma  followed  by  40  z  levels  extending  down  to 
5500  m.  The  initial  and  boundary  conditions  are  taken 
from  a  global  NCOM  model  with  Vs0  resolution.  Atmo¬ 
spheric  forcing  fields  (wind  stress,  atmospheric  pressure, 
and  heat  flux)  are  taken  from  the  Navy  Operational 
Global  Atmospheric  Prediction  System  (NO GAPS)  with 
0.5°  horizontal  resolution  (Rosmond  et  al.  2002). 

The  assimilation  system  used  here,  NCOM-4DVAR, 
is  dynamically  consistent  with  the  physics  and  numerics 
of  NCOM.  NCOM-4DVAR  uses  the  indirect  repre¬ 
senter  method  developed  by  Egbert  et  al.  (1994),  Chua 
and  Bennett  (2001),  and  Bennett  (2002).  See  Ngodock 
(2005)  and  Ngodock  and  Carrier  (2013)  for  more  details. 
The  velocity  and  mass  field  corrections  to  the  forward 
model  are  connected  by  a  barotropic  model  with  line¬ 
arized  covariance  matrix  operators. 

Here,  NCOM-4DVAR  is  run  in  a  weak-constraint 
mode,  so  that  possible  errors  in  the  model  dynamics  are 
included  in  the  assimilation,  in  addition  to  initial  con¬ 
dition  errors  (Ngodock  and  Carrier  2014).  The  initial 
errors  at  the  surface  are  set  to  2°C  for  temperature, 
0.5 psu  for  salinity,  0.2 ms-1  for  u  and  v  velocity,  and 
0.1m  for  SSH.  These  errors  are  attenuated  with  in¬ 
creasing  depth  (Yu  et  al.  2012).  This  specification  is 
appropriate  because  maximum  errors  are  expected  at 
the  surface  due  to  uncertainties  in  atmospheric  forcing. 
The  error  magnitudes  represent  10%  of  the  average 
surface  forcing  fields  over  the  assimilation  window.  This 
choice  is  consistent  with  that  of  Carrier  et  al.  (2014).  The 
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choice  of  10%  is  deemed  acceptable  as  ocean  surface 
currents  are  strongly  influenced  by  surface  wind  stress. 
The  horizontal  attenuation  length  scales  of  both  the 
initial  and  model  errors  range  from  30  to  40  km,  corre¬ 
sponding  to  the  Rossby  radius  of  deformation  in  the 
GOM.  While  the  background  error  scales  are  static 
through  time,  the  flow-dependent,  time-evolving  scales 
are  applied  by  the  dynamics  from  the  tangent  linear  and 
adjoint  models  of  NCOM. 

Each  observation  is  also  assigned  an  error  value. 
Observation  errors  include  both  estimated  instrument 
error  and  representative  error.  Representative  errors 
account  for  processes  unresolved  by  the  model.  A 
comparison  between  our  experiments  and  a  higher- 
resolution  (3  km)  reanalysis  (named  RT3km_v2  in 
Jacobs  et  al.  2014b)  shows  that  the  variance  of  the  dif¬ 
ference  of  surface  velocities  is  approximately  0.03  ms-1 
and  is  subsumed  into  the  larger  observational  error. 
Observation  errors  are  assumed  to  be  constant  in  time, 
with  values  of  0.2°C  for  temperature,  0.1  psu  for  salinity, 
and  0.05 ms-1  for  velocity.  These  temperature  and 
salinity  errors  are  representative  of  values  used  by  the 
operational  data  processing  system  used  by  the  U.S. 
Navy.  As  noted  by  Carrier  et  al.  (2014),  the  assumed 
velocity  observation  error  results  from  drifter  position 
uncertainties  (10  m)  and  the  time  interval  between 
positions  (15  min).  The  velocity  observation  error 
(0.02 ms-1)  was  conservatively  estimated  as  the  worst 
case,  with  maximum  error  in  both  drifter  positions 
totaling  20  m. 

3.  Observations  and  data  processing 

During  the  Grand  Lagrangian  Deployment  (GLAD), 
the  Consortium  for  Advanced  Research  on  Transport 
of  Hydrocarbons  in  the  Environment  (CARTHE) 
deployed  roughly  300  CODE-type  surface  drifters 
(drogued  at  1  m)  in  the  northern  GOM  at  the  end  of  July 
2012  (Poje  et  al.  2014).  These  drifters  were  tracked  using 
commercially  available  GPS  units,  and  reported  their 
positions  via  satellite  roughly  every  5  min. 

For  this  analysis,  GLAD  drifter  positions  were  as¬ 
sumed  to  be  accurate  to  within  10  m.  The  GPS  manu¬ 
facturer  specifies  an  accuracy  of  6.4  m.  A  dock-side  test 
verified  that  98%  of  reported  positions  were  within  10  m 
of  the  mean  measurement  for  each  of  10  drifters  over 
a  12-h  period.  Initial  position  uncertainty  of  this  mag¬ 
nitude  leads  to  endpoint  uncertainty  for  model  trajec¬ 
tories  after  96  h  on  the  order  of  200  m.  Thus,  both 
observed  and  modeled  position  uncertainties  are  negli¬ 
gible  relative  to  the  model  errors  discussed  below. 

A  drifter  position  record  was  terminated  at  the  point 
where  the  signal  was  lost  for  more  than  24  h,  if  the  drifter 


was  known  to  be  picked  up  by  a  boat,  or  if  it  traveled 
more  than  80  km  in  a  12-h  period  (implying  a  mean 
speed  of  1.85 ms-1  over  12 h). 

Infrequent  large  errors  were  found  in  the  GLAD 
drifter  position  records,  typically  due  to  poor  GPS  sat¬ 
ellite  reception.  These  outliers  were  identified  and  re¬ 
moved  by  a  three-step  procedure:  1)  positions  that 
implied  instantaneous  drifter  speeds  greater  than 
3 ms-1  were  deleted,  2)  positions  that  implied  a  drifter 
track  that  rotated  through  more  than  360  compass  de¬ 
grees  within  3  hours  were  deleted,  and  3)  each  remaining 
position  was  estimated  with  a  local  spline  interpolation 
using  past  and  future — but  not  the  current — positions 
from  the  same  drifter.  Measured  positions  that  deviated 
more  than  100  m  from  the  corresponding  spline  estimate 
were  deleted. 

When  all  outliers  were  removed,  record  gaps  were 
filled  by  spline  interpolation  of  the  remaining  valid  po¬ 
sitions  to  uniform  5-min  intervals.  Estimates  of  u  (east- 
west)  and  v  (north-south)  velocities  were  computed 
from  these  5-min  records  using  second-order,  centered 
finite  differences. 

To  reduce  noise,  the  5-min  position  and  velocity  re¬ 
cords  were  filtered  using  a  fourth-order  Butterworth 
low-pass  filter  with  a  1-h  period  cutoff.  Finally,  the  low- 
pass-filtered  records  were  subsampled  at  uniform 
1 5-min  intervals  beginning  on  whole  hours.  The  GLAD 
trajectory  dataset  used  here  is  publicly  available  (https:// 
data.gulfresearchinitiative.org/data/Rl.xl34.073:0004/). 
Figure  1  shows  a  time  series  of  the  number  of  GLAD 
trajectories  available  for  each  12-h  window  during  the 
study  period. 

Additional  assimilated  observations  include  GOES- 
East  sea  surface  temperature  (SST),  ARGO  profiling 
floats  (Roemmich  et  al.  2001),  and  expendable  bathy¬ 
thermographs  (XBT).  These  represent  typical  observa¬ 
tion  types  and  distributions  for  the  region  of  interest. 
The  quality  control  and  preparation  of  these  observa¬ 
tions  is  performed  by  the  NCODA  system  (Cummings 
2005).  Note  that  the  vast  majority  of  these  observations 
are  SST  data. 

The  time  period  for  the  analysis  presented  here  is  1 
August-30  September  2012.  During  this  period,  a  large 
Loop  Current  eddy  (evident  in  the  23  August  SSH  field 
shown  in  Fig.  2)  was  present  in  the  mid-GOM.  Figure  3 
shows  that  the  GLAD  trajectories  covered  much  of  the 
northeastern  GOM  during  this  period,  with  several 
preferred  pathways  evident.  Velocities  derived  from 
these  drifter  trajectories  help  to  constrain  the  mesoscale 
features  of  the  model  forecasts.  An  important  event 
during  the  study  period  was  the  passage  of  Hurricane 
Isaac,  a  category  1  hurricane,  during  26-30  August  2012. 
The  storm  was  responsible  for  the  sharp  drop  in  the 
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Fig.  1.  Number  of  GLAD  drifters  during  the  study  period. 


number  of  reporting  drifters  seen  in  Fig.  1.  It  also  spread 
the  remaining  drifters  over  a  much  larger  area  of  the 
GOM.  While  not  quantifiable,  it  is  possible  that  the 
loss  in  numbers  was  compensated  for  by  better  spatial 
coverage. 


4.  Methods 

In  this  section,  two  NCOM-4DVAR  model  forecast 
experiments  with  assimilation  are  compared  with  a  free- 
running  NCOM  model  without  any  assimilation.  The 
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Fig.  2.  Composite  GOM  SSH  field  showing  the  position  of  a  Loop  Current  ring  on  23  Aug  2012.  [Courtesy 
of  Colorado  Center  for  Astrodynamics  Research  (CCAR).] 
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Fig.  3.  GLAD  drifter  trajectories  for  the  period  1  Aug-30  Sep  2012. 


free  running  experiment  is  referred  to  as  the  “FREE” 
run.  The  first  assimilative  forecast  uses  only  temperature 
and  salinity  observations  and  it  is  referred  to  as  the 
“T/S”  run.  The  second  assimilative  forecast  uses 
temperature,  salinity,  and  velocities  from  the  GLAD 
drifters  and  it  is  referred  to  as  the  “ALL”  run.  The  two 
assimilative  forecasts  were  run  as  described  in  Carrier 
et  al.  (2014).  Here  “assimilative  forecasts”  refers  to 
those  forecasts  initialized  from  the  assimilation 
analysis. 

The  T/S  and  ALL  forecasts  use  4-day  assimilation 
windows  throughout  the  study  period.  After  each  4-day 
assimilation  cycle,  4-day  forecasts  are  produced  by 
the  forward  model  using  the  corrected  end  state  of 
the  previous  assimilation  cycle  as  the  initial  condi¬ 
tion.  The  4-day  forecasts  are  later  used  as  the  back¬ 
ground  for  the  next  4-day  assimilation  window.  As 
described  by  Carrier  et  al.  (2014),  data  injection  every 
hour  during  the  assimilation  period  produced  the  best 
fit  to  the  observations.  All  data  each  hour  was  spa¬ 
tially  subsampled  to  match  the  spatial  correlation 
length.  The  relatively  long  forecast  window  of  4  days 
was  chosen  to  achieve  a  reasonable  level  of  in¬ 
dependence  between  assimilated  and  evaluation 
data. 

Simulated  trajectories  are  computed  for  each  NCOM 
forecast  using  simple  advection.  No  attempt  is  made  to 
account  for  subgrid  scale  or  diffusion  effects.  Simulated 
drifters  are  advected  with  hourly  model  surface  veloci¬ 
ties  (linearly  interpolated  in  space  and  time)  using  an 
explicit,  adaptive  time-step  fourth-order  Runge-Kutta 
scheme.  The  initial  locations  of  the  simulated  drifters 


are  set  as  the  locations  of  the  GLAD  drifters  at  the  start 
time  of  each  4-day  forecast. 

Two  metrics  are  used  to  compare  GLAD  trajectories 
with  simulated  trajectories  from  all  three  NCOM  fore¬ 
casts.  The  first  metric  is  separation  distance ,  computed  as 
the  distance  (in  kilometers)  between  an  observed  and 
simulated  drifter  at  a  specified  time.  This  metric  is  fre¬ 
quently  used  to  assess  trajectory  forecast  skill  (e.g.,  see 
Huntley  et  al.  2011).  The  second  metric  is  angular  dif¬ 
ference ,  used  to  investigate  possible  direction  biases  in 
model  forecasts.  A  drifter  angle  is  defined  as  the  compass 
angle  of  the  line  drawn  between  the  drifter  launch  posi¬ 
tion  and  its  position  at  a  specified  time.  The  angular  dif¬ 
ference  is  defined  as  the  difference  (in  degrees)  between 
observed  and  simulated  drifter  angles.  These  metrics  are 
computed  using  all  available  GLAD  trajectories. 

5.  Results 

Figure  4  shows  a  time  series  of  mean  separation  dis¬ 
tance  (computed  using  all  available  GLAD  trajectories) 
24  h  into  the  4-day  forecast  for  all  three  NCOM  runs 
over  the  14  assimilation  intervals  in  the  study  period. 
Over  the  entire  study  period,  the  ALL  run  produces  the 
smallest  separation  distances,  with  a  mean  value  of  ap¬ 
proximately  10  km.  Separation  distance  for  the  ALL  run 
also  decreases  with  time  as  the  model  is  repeatedly 
corrected  with  the  observed  velocities.  Separation  dis¬ 
tances  after  48,  72,  and  96  h  (not  shown)  show  similar 
behavior  over  time,  with  larger  magnitudes.  It  is  also 
interesting  to  note  that  the  assimilation  of  just  temper¬ 
ature  and  salinity  does  not  outperform  the  FREE  run  in 
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Fig.  4.  Mean  separation  distance  (km)  after  24  h  between  the  GLAD  and  simulated  drifters  for 
each  assimilation  cycle.  The  average  is  taken  over  all  the  available  drifters  at  each  time. 


predicting  velocity  trajectories  (especially  toward  the 
later  part  of  the  experiment).  This  decrease  in  forecast 
skill,  due  to  the  low  predictability  of  the  ocean  response 
to  hurricanes,  coupled  with  the  ageostrophic  conditions 
that  persist  in  the  wake  of  a  passing  hurricane,  make  the 
constraining  of  surface  velocities  with  temperature  and 
salinity  observations  difficult. 

Separation  distances  for  the  FREE  run  and  T/S 
forecasts  peak  on  21  August,  after  the  fourth  assimila¬ 
tion  cycle  (Fig.  4).  For  the  ALL  forecast,  which  assimi¬ 
lates  GLAD  velocities,  this  peak  is  reduced  markedly. 
The  peak  on  21  August  is  due  to  forecast  errors  (most 
noticeable  in  the  T/S  and  FREE  runs)  in  the  position 
and  shape  of  the  Loop  Current  ring.  Figures  5a  and  5b 
show  the  mean  FREE  and  T/S  forecast  SSH  field  com¬ 
puted  over  the  period  21-25  August  2012,  as  well  as 
some  example  GLAD  and  model  trajectories.  There  is 
an  obvious  mismatch  between  the  GLAD  trajectories 
and  the  mean  FREE  and  T/S  forecast  SSH.  Figure  5c 
shows  the  mean  ALL  forecast  SSH  field  computed  over 
the  same  time  period.  The  ALL  forecast  produces 
a  better  realization  of  the  Loop  Current  ring  yielding 
much  better  agreement  between  the  GLAD  and  simu¬ 
lated  drifters.  For  the  ALL  forecast  (Fig.  5c)  both  the 
observed  and  simulated  trajectories  follow  the  ring 
boundary,  consistent  with  geostrophic  flow.  The  Loop 
Current  ring  shape  in  Fig.  5c  agrees  very  well  with  the 
composite  altimetry  map  for  23  August  shown  in  Fig.  2. 
The  assimilation  of  the  velocity  data  improves  the  SSH 
field  due  to  the  dynamical  connection  between  the  mass 
and  velocity  fields  provided  by  the  barotropic  model. 
The  extent  of  this  Eulerian  improvement  is  explored 
more  fully  in  a  companion  paper  (Carrier  et  al.  2014). 

Time  series  of  separation  distance,  computed  hourly 
as  the  average  over  all  available  drifters,  for  the  96-h 


forecast  of  each  assimilation  cycle  are  shown  for  the 
FREE  forecasts  (Fig.  6a),  the  T/S  forecasts  (Fig.  6b),  and 
the  ALL  forecasts  (Fig.  6c).  In  these  figures,  one  curve  is 
shown  for  each  of  the  14  assimilation  cycles  during  the 
study  period.  The  ALL  forecasts  (Fig.  6c)  consistently 
produce  better  trajectory  forecasts,  with  much  lower 
separation  distances  over  all  cycles  when  compared  to 
the  FREE  (Fig.  6a)  and  T/S  forecasts  (Fig.  6b).  When 
the  14  curves  in  Figs.  6a-c  are  averaged,  linear  fits  to  the 
curve  show  a  12.5  km  day-1  increase  in  separation  dis¬ 
tance  for  the  ALL  forecasts  compared  to  19.8  km  for  the 
T/S  forecasts  and  an  18.6  km  day-1  for  the  FREE  fore¬ 
casts.  Figure  6  also  shows  the  increased  spread  of  the 
separation  distance  enveloped  for  the  FREE  and  T/S 
forecasts  when  compared  to  the  ALL  forecasts. 

Time  series  of  the  distributions  of  angular  differences 
are  best  presented  as  two-dimensional  histograms.  His¬ 
tograms  of  angular  differences  after  24-h  forecasts  are 
shown  for  the  FREE  forecasts  (Fig.  7a),  T/S  forecasts 
(Fig.  7b),  and  the  ALL  forecasts  (Fig.  1C).  For  these 
histograms,  angular  differences  were  binned  in  15°  in¬ 
tervals.  A  comparison  among  Figs.  7a-c  shows  that  the 
ALL  forecasts  produce  smaller  angular  differences,  with 
most  angular  differences  ranging  from  approximately 
-60°  to  60°.  The  FREE  and  T/S  experiments  (Figs.  7a 
and  7b)  have  more  drifters  with  larger  angular  differ¬ 
ences.  The  histograms  in  Fig.  7  do  not  show  any  obvious 
preferred  directional  mismatch,  with  most  angular  dif¬ 
ferences  ranging  from  approximately  —45°  to  45°. 

To  assess  the  potential  spatial  range  of  our  model 
corrections  for  improving  Lagrangian  predictability,  a 
data-denial  experiment  was  performed,  where  approxi¬ 
mately  50%  of  the  original  drifter  record  was  left  out  for 
validation.  Figure  8  shows  the  trajectories  of  the  in¬ 
cluded  drifters  as  well  as  a  50-km-radius  circle.  Any 
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Fig.  5.  Mean  SSH  (m)  for  the  (a)  FREE  experiment,  (b)  T/S  experiment,  and  (c)  ALL 
experiment  for  0000  UTC  21  Aug-0000  UTC  25  Aug  2012.  Black  crosses  indicate  observed 
positions  at  0000  UTC  21  Aug  2012  for  six  GLAD  drifters.  Observed  (green)  and  simulated 
(purple)  trajectories  using  (a)  FREE,  (b)  T/S,  and  (c)  ALL  forecast  velocities  are  also  shown. 


drifter  that  entered  that  circle  during  the  2-month  study 
period  was  removed  from  the  assimilation.  This  location 
was  chosen  to  make  this  test  as  rigorous  as  possible  by 
removing  a  large  number  of  drifters  that  observe  the 


mesoscale  eddy  that  is  a  main  circulation  feature  in  this 
domain. 

Figure  9a  shows  the  separation  distance  time  series 
similar  to  Fig.  4  except  that  the  separation  distance  is 
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Fig.  6.  Mean  separation  distance  as  a  function  of  forecast  length 
for  each  of  the  14  assimilation  cycles  during  the  two-month  study 
period  (gray  lines)  for  the  (a)  FREE  experiment,  (b)  T/S  experi¬ 
ment,  and  (c)  ALL  experiment,  with  the  linear  fit  to  all  the  cycles 
(black  line).  The  average  is  taken  over  all  the  available  drifters  for 
each  time  step. 

calculated  only  at  the  locations  of  the  removed  drifters 
and  the  data-denial  experiment  is  included  as  well.  The 
data-denial  experiment  produces  significantly  lower 
separation  distances  than  the  FREE  and  T/S  runs  until 
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Fig.  7.  Two-dimensional  histograms  showing  the  percentage  of 
drifters  with  a  given  mean  angular  model-observation  mismatch 
after  24  h  for  each  assimilation  cycle,  (a)  FREE  experiment,  (b)  T/S 
experiment,  and  (c)  ALL  experiment. 

the  beginning  of  September.  Thereafter  the  separation 
distances  for  the  data-denial  experiment  are  of  the  same 
order  as  the  FREE  and  T/S  runs,  with  the  ALL  experi¬ 
ment  producing  the  lowest  distances.  As  mentioned 
earlier  the  removed  drifters  were  in  the  region  of  the 
large  eddy  whose  misplacement  causes  the  large  peak  in 
the  separation  distances  of  the  FREE  and  T/S  runs. 
Even  though  the  majority  of  the  remaining  drifters  did 
not  directly  sample  this  eddy,  the  data-denial  experi¬ 
ment  still  successfully  reduces  the  separation  distances 
because  sufficient  drifters  remained  to  constrain  this 
eddy  feature. 
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Fig.  8.  GLAD  drifter  trajectories  assimilated  for  the  period  1  Aug-30  Sep  2012  for  the  data- 
denial  experiment.  All  drifters  that  ever  entered  the  black  circle  were  removed  from  the 
original  record. 


The  improvement  of  the  Lagrangian  predictability  for 
the  data-denial  experiment  over  the  FREE  and  T/S  runs 
disappears  in  September  (Fig.  9a).  This  result  can  be 
explained  by  the  fact  that  in  this  month  the  distance 
between  unassimilated  and  assimilated  drifters  signifi¬ 
cantly  increased.  Figure  9b  shows  the  time  series  of  the 
mean  over  all  unassimilated  drifters  of  the  distance  to 
the  nearest  assimilated  drifter  for  each  time  step.  The 
large  peak  near  the  end  of  August  coincides  with  the 
failure  of  a  single  assimilated  drifter  on  27  August,  which 
had  been  constraining  a  large  number  of  unassimilated 
drifters.  The  hurricane  passage  around  the  same  time 
also  spread  the  drifters  out  spatially.  Thus,  the  period  of 
improved  predictability  (9-27  August)  is  marked  by  an 
average  distance  to  an  assimilated  drifter  of  roughly  15- 
40  km,  whereas  during  the  following  period  (28  August- 
29  September)  the  average  distance  to  an  assimilated 
drifter  is  much  larger,  around  60-80  km. 

These  results  suggest  that  a  critical  distance  above 
which  an  assimilated  drifter  will  not  sufficiently  con¬ 
strain  an  unassimilated  one  is  around  40-60  km.  It 
should  be  cautioned,  however,  that  this  critical  distance 
is  likely  to  depend  strongly  on  the  flow  regime  and  the 
dominant  scales  of  circulation  features  in  the  area  of 


interest.  Moreover,  it  is  likely  to  be  nonisotropic. 
Clearly,  more  analysis  is  needed  for  a  more  compre¬ 
hensive  assessment  of  the  spatial  scales  of  a  drifter’s 
impact  on  nearby  trajectories,  which  we  hope  to  report 
on  in  follow-on  papers. 

6.  Summary 

A  previous  study  (Carrier  et  al.  2014)  showed  that 
assimilating  velocity  data  inferred  from  drifters  using 
4DVAR  can  markedly  improve  model  temperature, 
salinity,  SSH,  and  velocity  Eulerian  forecast  skill.  Here, 
analysis  of  the  same  forecasts  described  by  Carrier  et  al. 
(2014)  shows  that  assimilating  these  drifter  velocities 
also  markedly  improves  trajectory  forecast  skill  over 
a  24-h  period  (see  Fig.  4).  Averages  over  the  fourteen 
4-day  assimilation  cycles  studied  here  show  that  assim¬ 
ilating  drifter  velocities  reduces  trajectory  forecast  error 
growth  rates  by  7.3  km  day-1,  a  substantial  improvement 
(see  Fig.  6).  Trajectory  direction  differences  also  de¬ 
crease  markedly  when  drifter  velocities  are  assimilated. 
The  data-denial  experiment  is  a  first  step  toward  de¬ 
termining  the  spatial  extent  of  the  improvements  to 
Lagrangian  predictability.  The  results  show  that  the 
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Fig.  9.  (a)  Mean  separation  distance  (km)  after  24  h  between  the  GLAD  and  simulated  drifters  for  each 
assimilation  cycle.  The  average  is  taken  only  over  the  drifters  removed  in  the  data-denial  experiment  (cf.  Fig.  4). 
(b)  The  minimum  distance  of  each  unassimilated  drifter  to  the  nearest  assimilated  drifter,  averaged  over  all  the 
unassimilated  drifters  at  each  hourly  time  step. 


Lagrangian  predictability  significantly  falls  off  away 
from  the  assimilated  drifters  for  distances  greater  than 
about  40  km.  Note  that  the  correlation  scale  selected  for 
these  Gulf  of  Mexico  experiments  is  also  40  km. 

Here,  the  dramatic  improvements  in  mean  trajec¬ 
tory  forecast  skill  were  due  to  increased  accuracy  of 
the  mesoscale  forecast  position  and  shape  of  a  large 
Loop  Current  ring  in  the  mid-GOM  (see  Fig.  5).  The 
added  constraints  provided  by  assimilated  drifter  ve¬ 
locities  not  only  improved  forecast  trajectories,  but 
also  provided  important  corrections  to  the  forecast 
SSH  field. 

The  tracking  and  prediction  of  oil  transport  at  the 
ocean  surface  remains  a  challenging  problem  for  ocean¬ 
ographers.  The  assimilation  of  these  high-resolution 
drifter  velocities  is  shown  to  provide  an  important  cor¬ 
rection  to  the  model  mesoscale  forecast  that  leads  to 


a  substantial  improvement  in  trajectory  forecast  skill.  This 
result  underscores  the  importance  of  drifter  observations, 
and  their  ability  to  substantially  improve  model  trajectory 
forecasts  in  response  to  future  oil  spills. 
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